#include <simgear/debug/logstream.hxx>

#include "tg_polygon_set.hxx"

#include <simgear/misc/sg_path.hxx> // for file i/o

// we are loading polygonal data from untrusted sources
// high probability this will crash CGAL if we just load 
// the points.  previous terragear would attempt to clean
// the input using duplicate point detection, degenerate
// edge detection, etc.
// 
// Instead, we'll generate an arrangement from each ring
// the first ring is considered the boundary, and all of 
// the rest are considered holes.
// this should never fail.                   _
//           _______________________________|_|
//          /                            ___|
//         /    _______                 / 
//        /    |       |               /
//        |    |     __|_             /
//        |    |____|__| |     ______/____
//        |         |____|     |____/     |
//        |                   /|          |
//        \                  / |__________|
//         \________________/
//
// NOTE: using this method, the above shapefile will result in a
// single polygon with holes.
// 1) the outer boundary is not simple - it self intersects in the 
//    top right corner.  We generate the outer boundary as the 
//    union of all faces generated by the first ring
// 2) the three remaining rings are unioned together as holes
//    a boolean difference is performed to make them holes.
//
// NOTE: the first two intersecting holes become a single hole.
//       the third ring decreases the boundary of the polygon
//
// the final result is two polygons_with_holes.
// the frst is a poly with a single hole.
// the second is the degenerate piece in the top right.
tgPolygonSet::tgPolygonSet( OGRPolygon* poGeometry, const tgPolygonSetMeta& metaInfo ) : meta(metaInfo)
{
    std::vector<cgalPoly_Polygon>	boundaries;
    std::vector<cgalPoly_Polygon>	holes;
    cgalPoly_PolygonSet             holesUnion;
    std::vector<cgalPoly_Point>     nodes;

    SG_LOG( SG_GENERAL, SG_DEBUG, "Geometry has : " <<  poGeometry->getNumInteriorRings() << " rings" );
    
    // create PolygonSet from the outer ring
    OGRLinearRing const *ring = poGeometry->getExteriorRing();
    nodes.clear();
    for (int i = 0; i < ring->getNumPoints(); i++) {
        nodes.push_back( cgalPoly_Point(ring->getX(i), ring->getY(i)) );
    }
    facesFromUntrustedNodes( nodes, boundaries, holes );

    SG_LOG( SG_GENERAL, SG_DEBUG, "Outer boundary complete : boundaries: " <<  boundaries.size() << " holes: " << holes.size() );

    // then a PolygonSet from each interior ring
    for ( int i = 0 ; i < poGeometry->getNumInteriorRings(); i++ ) {
        ring = poGeometry->getInteriorRing( i );
        nodes.clear();
        for (int j = 0; j < ring->getNumPoints(); j++) {
            nodes.push_back( cgalPoly_Point(ring->getX(j), ring->getY(j)) );
        }
        facesFromUntrustedNodes( nodes, holes, boundaries );
        
        SG_LOG( SG_GENERAL, SG_DEBUG, "hole " << i << " of " << poGeometry->getNumInteriorRings() << " complete : boundaries: " <<  boundaries.size() << " holes: " << holes.size() );        
    }
    
    // join all the boundaries
    ps.join( boundaries.begin(), boundaries.end() );
    
    // join all the holes
    holesUnion.join( holes.begin(), holes.end() );
    
    // perform difference
    ps.difference( holesUnion );

    if ( ps.is_empty() ) {
        SG_LOG( SG_GENERAL, SG_ALERT, "tgPolygonSet::tgPolygonSet DIFF between boundaries and holes is empty" );
    }
}

tgPolygonSet::tgPolygonSet( OGRFeature* poFeature, OGRPolygon* poGeometry )
{
    std::vector<cgalPoly_Polygon>	boundaries;
    std::vector<cgalPoly_Polygon>	holes;
    cgalPoly_PolygonSet             holesUnion;
    std::vector<cgalPoly_Point>     nodes;

    // generate texture info from feature
    meta.getFeatureFields( poFeature );
    
    SG_LOG( SG_GENERAL, SG_DEBUG, "got material: " << meta.getMaterial() );
    
    // create PolygonSet from the outer ring
    OGRLinearRing const *ring = poGeometry->getExteriorRing();
    nodes.clear();
    for (int i = 0; i < ring->getNumPoints(); i++) {
        nodes.push_back( cgalPoly_Point(ring->getX(i), ring->getY(i)) );
    }
    facesFromUntrustedNodes( nodes, boundaries, holes );

    // then a PolygonSet from each interior ring
    for ( int i = 0 ; i < poGeometry->getNumInteriorRings(); i++ ) {
        ring = poGeometry->getInteriorRing( i );
        nodes.clear();
        for (int j = 0; j < ring->getNumPoints(); j++) {
            nodes.push_back( cgalPoly_Point(ring->getX(j), ring->getY(j)) );
        }
        facesFromUntrustedNodes( nodes, holes, boundaries );
    }

    // join all the boundaries
    ps.join( boundaries.begin(), boundaries.end() );

    // join all the holes
    holesUnion.join( holes.begin(), holes.end() );

    // perform difference
    ps.difference( holesUnion );
}

GDALDataset* tgPolygonSet::openDatasource( const char* datasource_name )
{
    GDALDataset*    poDS = NULL;
    GDALDriver*     poDriver = NULL;
    const char*     format_name = "ESRI Shapefile";
    
    SG_LOG( SG_GENERAL, SG_DEBUG, "Open Datasource: " << datasource_name );
    
    GDALAllRegister();
    
    poDriver = GetGDALDriverManager()->GetDriverByName( format_name );
    if ( poDriver ) {    
        poDS = poDriver->Create( datasource_name, 0, 0, 0, GDT_Unknown, NULL );
    }
    
    return poDS;
}

OGRLayer* tgPolygonSet::openLayer( GDALDataset* poDS, OGRwkbGeometryType lt, PolygonSetLayerFields lf, const char* layer_name )
{
    OGRLayer*           poLayer = NULL;
 
    if ( !strlen( layer_name )) {
        SG_LOG(SG_GENERAL, SG_ALERT, "tgPolygonSet::toShapefile: layer name is NULL" );
        exit(0);
    }
    
    poLayer = poDS->GetLayerByName( layer_name );    
    if ( !poLayer ) {
        SG_LOG(SG_GENERAL, SG_DEBUG, "tgPolygonSet::toShapefile: layer " << layer_name << " doesn't exist - create" );

        OGRSpatialReference srs;
        srs.SetWellKnownGeogCS("WGS84");
        
        poLayer = poDS->CreateLayer( layer_name, &srs, lt, NULL );

        OGRFieldDefn descriptionField( "tg_desc", OFTString );
        descriptionField.SetWidth( 128 );
        if( poLayer->CreateField( &descriptionField ) != OGRERR_NONE ) {
            SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_desc' failed" );
        }
        
        OGRFieldDefn idField( "tg_id", OFTInteger );
        if( poLayer->CreateField( &idField ) != OGRERR_NONE ) {
            SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_id' failed" );
        }

        if ( lf == LF_ALL ) {
            OGRFieldDefn metaField( "tg_meta", OFTInteger );
            if( poLayer->CreateField( &metaField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'meta' failed" );
            }

            OGRFieldDefn fidField( "OGC_FID", OFTInteger );
            if( poLayer->CreateField( &fidField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'OGC_FID' failed" );
            }

            OGRFieldDefn flagsField( "tg_flags", OFTInteger );
            if( poLayer->CreateField( &flagsField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'flags' failed" );
            }
            
            OGRFieldDefn materialField( "tg_mat", OFTString );
            materialField.SetWidth( 32 );
            if( poLayer->CreateField( &materialField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_material' failed" );
            }
            
            OGRFieldDefn texMethodField( "tg_texmeth", OFTInteger );
            if( poLayer->CreateField( &texMethodField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tex_method' failed" );
            }
            
            OGRFieldDefn texRefLonField( "tg_reflon", OFTReal );
            texRefLonField.SetWidth( 24 );
            texRefLonField.SetPrecision( 3 );        
            if( poLayer->CreateField( &texRefLonField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tp_ref_lon' failed" );
            }
            
            OGRFieldDefn texRefLatField( "tg_reflat", OFTReal );
            if( poLayer->CreateField( &texRefLatField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tp_ref_lat' failed" );
            }
            
            OGRFieldDefn texHeadingField( "tg_heading", OFTReal );
            if( poLayer->CreateField( &texHeadingField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tp_heading' failed" );
            }
            
            OGRFieldDefn texWidthField( "tg_width", OFTReal );
            if( poLayer->CreateField( &texWidthField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tp_width' failed" );
            }
            
            OGRFieldDefn texLengthField( "tg_length", OFTReal );
            if( poLayer->CreateField( &texLengthField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tp_length' failed" );
            }
            
            OGRFieldDefn texMinUField( "tg_minu", OFTReal );
            if( poLayer->CreateField( &texMinUField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tp_minu' failed" );
            }
            
            OGRFieldDefn texMinVField( "tg_minv", OFTReal );
            if( poLayer->CreateField( &texMinVField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tp_minv' failed" );
            }
            
            OGRFieldDefn texMaxUField( "tg_maxu", OFTReal );
            if( poLayer->CreateField( &texMaxUField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tp_maxu' failed" );
            }
            
            OGRFieldDefn texMaxVField( "tg_maxv", OFTReal );
            if( poLayer->CreateField( &texMaxVField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tp_maxv' failed" );
            }
            
            OGRFieldDefn texMinClipUField( "tg_mincu", OFTReal );
            if( poLayer->CreateField( &texMinClipUField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tp_min_clipu' failed" );
            }
            
            OGRFieldDefn texMinClipVField( "tg_mincv", OFTReal );
            texMinClipVField.SetWidth( 24 );
            texMinClipVField.SetPrecision( 3 );        
            if( poLayer->CreateField( &texMinClipVField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tp_min_clipv' failed" );
            }
            
            OGRFieldDefn texMaxClipUField( "tg_maxcu", OFTReal );
            if( poLayer->CreateField( &texMaxClipUField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tp_max_clipu' failed" );
            }
            
            OGRFieldDefn texMaxClipVField( "tg_maxcv", OFTReal );
            if( poLayer->CreateField( &texMaxClipVField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tp_max_clipv' failed" );
            }
            
            OGRFieldDefn texCenterLatField( "tg_clat", OFTReal );
            if( poLayer->CreateField( &texCenterLatField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tg_tp_center_lat' failed" );
            }


            // surface metadata
            OGRFieldDefn texSurfaceMinLonField( "tgsrf_mnln", OFTReal );
            texSurfaceMinLonField.SetWidth( 24 );
            texSurfaceMinLonField.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceMinLonField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_mnln' failed" );
            }

            OGRFieldDefn texSurfaceMinLatField( "tgsrf_mnlt", OFTReal );
            texSurfaceMinLatField.SetWidth( 24 );
            texSurfaceMinLatField.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceMinLatField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_mnlt' failed" );
            }

            OGRFieldDefn texSurfaceMaxLonField( "tgsrf_mxln", OFTReal );
            texSurfaceMaxLonField.SetWidth( 24 );
            texSurfaceMaxLonField.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceMaxLonField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_mxln' failed" );
            }

            OGRFieldDefn texSurfaceMaxLatField( "tgsrf_mxlt", OFTReal );
            texSurfaceMaxLatField.SetWidth( 24 );
            texSurfaceMaxLatField.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceMaxLatField ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_mxlt' failed" );
            }
            
            OGRFieldDefn texSurfaceCoef00( "tgsrf_co00", OFTReal );
            texSurfaceCoef00.SetWidth( 24 );
            texSurfaceCoef00.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceCoef00 ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_co00' failed" );
            }

            OGRFieldDefn texSurfaceCoef01( "tgsrf_co01", OFTReal );
            texSurfaceCoef01.SetWidth( 24 );
            texSurfaceCoef01.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceCoef01 ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_co01' failed" );
            }

            OGRFieldDefn texSurfaceCoef02( "tgsrf_co02", OFTReal );
            texSurfaceCoef02.SetWidth( 24 );
            texSurfaceCoef02.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceCoef02 ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_co02' failed" );
            }

            OGRFieldDefn texSurfaceCoef03( "tgsrf_co03", OFTReal );
            texSurfaceCoef03.SetWidth( 24 );
            texSurfaceCoef03.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceCoef03 ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_co03' failed" );
            }

            OGRFieldDefn texSurfaceCoef04( "tgsrf_co04", OFTReal );
            texSurfaceCoef04.SetWidth( 24 );
            texSurfaceCoef04.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceCoef04 ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_co04' failed" );
            }

            OGRFieldDefn texSurfaceCoef05( "tgsrf_co05", OFTReal );        
            texSurfaceCoef05.SetWidth( 24 );
            texSurfaceCoef05.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceCoef05 ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_co05' failed" );
            }

            OGRFieldDefn texSurfaceCoef06( "tgsrf_co06", OFTReal );
            texSurfaceCoef06.SetWidth( 24 );
            texSurfaceCoef06.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceCoef06 ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_co06' failed" );
            }
            
            OGRFieldDefn texSurfaceCoef07( "tgsrf_co07", OFTReal );
            texSurfaceCoef07.SetWidth( 24 );
            texSurfaceCoef07.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceCoef07 ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_co07' failed" );
            }

            OGRFieldDefn texSurfaceCoef08( "tgsrf_co08", OFTReal );
            texSurfaceCoef08.SetWidth( 24 );
            texSurfaceCoef08.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceCoef08 ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_co08' failed" );
            }
            
            OGRFieldDefn texSurfaceCoef09( "tgsrf_co09", OFTReal );
            texSurfaceCoef09.SetWidth( 24 );
            texSurfaceCoef09.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceCoef09 ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_co09' failed" );
            }

            OGRFieldDefn texSurfaceCoef10( "tgsrf_co10", OFTReal );
            texSurfaceCoef10.SetWidth( 24 );
            texSurfaceCoef10.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceCoef10 ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_co10' failed" );
            }

            OGRFieldDefn texSurfaceCoef11( "tgsrf_co11", OFTReal );
            texSurfaceCoef11.SetWidth( 24 );
            texSurfaceCoef11.SetPrecision( 3 );        
            if( poLayer->CreateField( &texSurfaceCoef11 ) != OGRERR_NONE ) {
                SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'tgsrf_co11' failed" );
            }
        }
    } else {
        SG_LOG(SG_GENERAL, SG_DEBUG, "tgPolygonSet::toShapefile: layer " << layer_name << " already exists - open" );        
    }
   
    return poLayer;
}

//void tgPolygonSet::toShapefile( OGRLayer* layer, const char* description ) const
//{
//    
//}

void tgPolygonSet::toShapefile( const char* datasource, const char* layer ) const
{
    // Open datasource and layer
    GDALDataset* poDS = openDatasource( datasource );

    if ( poDS ) {
        OGRLayer* poLayer = openLayer( poDS, wkbPolygon25D, LF_ALL, layer );
        //OGRLayer* poLayer = openLayer( poDS, wkbLineString25D, LF_ALL, layer );
        
        if ( poLayer ) {
            toShapefile( poLayer, ps );
        }
    }
    
    // close datasource
    GDALClose( poDS );
}

void tgPolygonSet::toShapefile( OGRLayer* poLayer ) const
{
    toShapefile( poLayer, ps );
}

void tgPolygonSet::toShapefile( OGRLayer* poLayer, const cgalPoly_PolygonSet& polySet ) const
{
    std::list<cgalPoly_PolygonWithHoles>                 pwh_list;
    std::list<cgalPoly_PolygonWithHoles>::const_iterator it;

    polySet.polygons_with_holes( std::back_inserter(pwh_list) );
    SG_LOG(SG_GENERAL, SG_DEBUG, "tgPolygonSet::toShapefile: got " << pwh_list.size() << " polys with holes ");
    
    // save each poly with holes to the layer
    for (it = pwh_list.begin(); it != pwh_list.end(); ++it) {
        cgalPoly_PolygonWithHoles pwh = (*it);

        toShapefile( poLayer, pwh );
    }
}

void tgPolygonSet::toShapefile( OGRLayer* poLayer, const cgalPoly_PolygonWithHoles& pwh ) const
{
    OGRPolygon    polygon;
    OGRPoint      point;
    OGRLinearRing ring;
    
    // in CGAL, the outer boundary is counter clockwise - in GDAL, it's expected to be clockwise
    cgalPoly_Polygon poly;
    cgalPoly_Polygon::Vertex_iterator it;

    poly = pwh.outer_boundary();
    //poly.reverse_orientation();
    for ( it = poly.vertices_begin(); it != poly.vertices_end(); it++ ) {
        point.setX( CGAL::to_double( (*it).x() ) );
        point.setY( CGAL::to_double( (*it).y() ) );
        point.setZ( 0.0 );
                
        ring.addPoint(&point);
    }
    ring.closeRings();
    polygon.addRing(&ring);

    // then write each hole
    cgalPoly_PolygonWithHoles::Hole_const_iterator hit;
    for (hit = pwh.holes_begin(); hit != pwh.holes_end(); ++hit) {
        OGRLinearRing hole;
        poly = (*hit);

        //poly.reverse_orientation();
        for ( it = poly.vertices_begin(); it != poly.vertices_end(); it++ ) {
            point.setX( CGAL::to_double( (*it).x() ) );
            point.setY( CGAL::to_double( (*it).y() ) );
            point.setZ( 0.0 );
                    
            hole.addPoint(&point);            
        }
        hole.closeRings();
        polygon.addRing(&hole);
    }

    OGRFeature* poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
    poFeature->SetGeometry(&polygon);
    meta.setFeatureFields( poFeature );
    
    if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE )
    {
        SG_LOG(SG_GENERAL, SG_ALERT, "Failed to create feature in shapefile");
    }
    OGRFeature::DestroyFeature(poFeature);    
}

void tgPolygonSet::toShapefile( OGRLayer* poLayer, const cgalPoly_Polygon& poly ) const
{
    OGRPolygon    polygon;
    OGRPoint      point;
    OGRLinearRing ring;
    
    // in CGAL, the outer boundary is counter clockwise - in GDAL, it's expected to be clockwise
    cgalPoly_Polygon::Vertex_iterator it;
    
    for ( it = poly.vertices_begin(); it != poly.vertices_end(); it++ ) {
        point.setX( CGAL::to_double( (*it).x() ) );
        point.setY( CGAL::to_double( (*it).y() ) );
        point.setZ( 0.0 );
        
        ring.addPoint(&point);
    }
    ring.closeRings();
    polygon.addRing(&ring);
    
    OGRFeature* poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
    poFeature->SetGeometry(&polygon);    
    meta.setFeatureFields( poFeature );
    
    if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE )
    {
        SG_LOG(SG_GENERAL, SG_ALERT, "Failed to create feature in shapefile");
    }
    OGRFeature::DestroyFeature(poFeature);    
}


void tgPolygonSet::toShapefile( OGRLayer* poLayer, const cgalPoly_Arrangement& arr ) const
{    
    cgalPoly_EdgeConstIterator eit;
    
    for ( eit = arr.edges_begin(); eit != arr.edges_end(); ++eit ) {        
        cgalPoly_Segment seg = eit->curve();

        OGRLinearRing ring;
        OGRPoint      point;

        point.setX( CGAL::to_double( seg.source().x() ) );
        point.setY( CGAL::to_double( seg.source().y() ) );
        point.setZ( 0 );
        ring.addPoint(&point);

        point.setX( CGAL::to_double( seg.target().x() ) );
        point.setY( CGAL::to_double( seg.target().y() ) );
        point.setZ( 0 );
        ring.addPoint(&point);

        OGRFeature* poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
        poFeature->SetGeometry(&ring);    
    
        if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE )
        {
            SG_LOG(SG_GENERAL, SG_ALERT, "Failed to create feature in shapefile");
        }
        OGRFeature::DestroyFeature(poFeature);    
    }    
}

void tgPolygonSet::fromShapefile( const OGRFeatureDefn* poFDefn, OGRCoordinateTransformation* poCT, OGRFeature* poFeature, tgPolygonSetList& polys )
{
    OGRGeometry *poGeometry = poFeature->GetGeometryRef();
    if (poGeometry == NULL) {
        SG_LOG( SG_GENERAL, SG_INFO, "Found feature without geometry!" );
        return;
    }
    
    OGRwkbGeometryType geoType = wkbFlatten(poGeometry->getGeometryType());
    if (geoType != wkbPolygon && geoType != wkbMultiPolygon) {
        SG_LOG( SG_GENERAL, SG_INFO, "Unknown feature" << geoType );
        return;
    }

    // now the geometry(s)
    poGeometry->transform( poCT );
            
    switch( geoType ) {
        case wkbPolygon:
        {
            tgPolygonSet poly( poFeature, (OGRPolygon *)poGeometry );
            polys.push_back(poly);
            break;
        }
        
        case wkbMultiPolygon:
        {
            OGRMultiPolygon* multipoly = (OGRMultiPolygon*)poGeometry;
            SG_LOG( SG_GENERAL, SG_DEBUG, "loaded multi poly with " << multipoly->getNumGeometries() << " polys" );
            
            for (int i=0;i<multipoly->getNumGeometries();i++) {
                tgPolygonSet poly( poFeature, (OGRPolygon *)multipoly->getGeometryRef(i) );
                polys.push_back(poly);
            }
            break;
        }
        
        default:
            break;
    }
    
    return;
}

void tgPolygonSet::processLayer(OGRLayer* poLayer, tgPolygonSetList& polys )
{
    OGRFeatureDefn*                 poFDefn = NULL;
    OGRCoordinateTransformation*    poCT = NULL;
    char*                           srsWkt;
    OGRSpatialReference*            oSourceSRS;
    OGRSpatialReference             oTargetSRS;
    OGRFeature*                     poFeature = NULL;
    std::string                     layername;
        
    /* determine the indices of the required columns */
    poFDefn = poLayer->GetLayerDefn();
    layername = poFDefn->GetName();
    
    /* setup a transformation to WGS84 */
    oSourceSRS = poLayer->GetSpatialRef();
    if (oSourceSRS == NULL) {
        SG_LOG( SG_GENERAL, SG_ALERT, "Layer " << layername << " has no defined spatial reference system" );
        exit( 1 );
    }
    
    oSourceSRS->exportToWkt(&srsWkt);
    SG_LOG( SG_GENERAL, SG_DEBUG, "Layer: " << layername << " spatial reference system: " << srsWkt );
    CPLFree(srsWkt);
    
    oTargetSRS.SetWellKnownGeogCS( "WGS84" );    
    poCT = OGRCreateCoordinateTransformation(oSourceSRS, &oTargetSRS);

    SG_LOG( SG_GENERAL, SG_DEBUG, "Layer " << layername << " has " << poLayer->GetFeatureCount() << " features " );
    
    // Generate the work queue for this layer
    while ( ( poFeature = poLayer->GetNextFeature()) != NULL )
    {
        fromShapefile( poFDefn, poCT, poFeature, polys );
        OGRFeature::DestroyFeature( poFeature );
    }
    
    OCTDestroyCoordinateTransformation ( poCT );
}

void tgPolygonSet::fromShapefile( const SGPath& p, tgPolygonSetList& polys )
{
    GDALDataset* poDS = NULL;
    OGRLayer*    poLayer = NULL;
        
    GDALAllRegister();
    
    poDS = (GDALDataset*)GDALOpenEx( p.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL );
    if( poDS == NULL )
    {
        SG_LOG( SG_GENERAL, SG_DEBUG, "Failed opening datasource " << p.c_str() );
        return;
    }
    
    SG_LOG( SG_GENERAL, SG_DEBUG, "Processing datasource " << p.c_str() << " with " << poDS->GetLayerCount() << " layers " );
    polys.clear();
    
    for (int i=0; i<poDS->GetLayerCount(); i++) {
        poLayer = poDS->GetLayer(i);
            
        assert(poLayer != NULL);
        processLayer(poLayer, polys );
    }
    
    GDALClose( poDS );    
    
    for ( unsigned int i=0; i<polys.size(); i++ ) {
        SG_LOG( SG_GENERAL, SG_DEBUG, "return poly " << i << " with material " << polys[i].getMeta().getMaterial() );
    }        
}
